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The quasi-particle energy spectrum of the Hartree-Fock-Bogoliubov (HFB) equations contains 
discrete bound states, resonances, and non-resonant continuum states. We study the structure of 
the unbound quasi-particle spectrum of weakly bound nuclei within several methods that do not 
rely on imposing scattering or outgoing boundary conditions. Various approximations are examined 
to estimate resonance widths. It is shown that the stabilization method works well for all HFB 
resonances except for very narrow ones. The Thomas-Fermi approximation to the non-resonant 
continuum has been shown to be very effective, especially for coordinate-space HFB calculations in 
large boxes that involve huge amounts of discretized quasi-particle continuum states. 
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I. INTRODUCTION 

A major challenge for theoretical nuclear structure re- 
search is the development of robust models and tech- 
niques aiming at the microscopic description of the nu- 
clear many-body problem and capable of extrapolating 
into unknown regions of the nuclear landscape. Since 
the coherent theoretical framework should describe both 
well-bound and drip-line nuclei, of particular importance 
is the treatment of the particle continuum in weakly 
bound nuclei and the development of microscopic reac- 
tion theory that is integrated with improved structure 
models [l|- 

For open-shell medium-mass and heavy nuclei, the the- 
oretical tool of choice is the nuclear Density Functional 
Theory (DFT) [2| . Its main ingredient is the energy den- 
sity functional that depends on proton and neutron den- 
sities and currents, as well as pairing densities describ- 
ing nuclear superconductivity "^. The HFB equations 
of nuclear DFTproperly take into account the scatter- 
ing continuum 1^-ul , and this is a welcome feature of the 
formalism. 

The quasi-particle energy spectrum of HFB consists 
of a finite number of bound quasi-particle states and a 
continuum of unbound states. The so-called deep-hole 
quasi-particle resonances are unique to the HFB theory 
[6j . They can be associated with single-particle Hartree- 
Fock (HF) states that are well bound in the absence of 
pairing correlations but acquire a finite particle width 
due to the continuum coupling induced by the particle- 
particle channel. As discussed in the literature 0, Q, 
those resonances cannot be properly described by BCS- 
like theories [S, 10] which yield particle and pairing den- 
sities which are not localized in space. 

The role played by the HFB continuum has been a 
subject of many works [3t3j Iii1'21]. The proper treat- 
ment of quasi-particle continuum is important not only 
for ground-state properties but also for the description of 



nuclear excitations, e.g., within a self-consistent QRPA 
approach ;22-26j. Within the real-energy HFB frame- 
work, the proper theoretical treatment of the HFB con- 
tinuum is fairly sophisticated since the scattering bound- 
ary conditions must be met. 

If the outgoing boundary conditions are imposed, 
the unbound HFB eigenstates have complex energies; 
their imaginary parts are related to the particle width. 
The complex-energy spherical HFB equations have been 
solved in Ref. Q within the Gamow HFB (GHFB) ap- 
proach, which shares many techniques with the Gamow 
shell model [13, [ll]. Within GHFB, quasi-particle reso- 
nance widths can be calculated with a high precision. Al- 
ternatively, diagonalizing the HFB matrix in the Gamow 
HF or Poschl-Teller-Ginocchio basis turned out to be an 
efficient way to account for the continuum effects [29] . 

In addition to the methods that employ correct asymp- 
totic boundary conditions for unbound states, the quasi- 
particle continuum can be approximately treated by 
means of a discretization method. The commonly used 
ap pro ach is to impose the box boundary conditions 
[3, ^^ [30l - l34| , in which wave functions are spanned by 
a basis of orthonormal functions defined on a lattice in 
coordinate space and enforced to be zero at box bound- 
aries. In this way, referred to as the C^ discretization, 
quasi-particle continuum of HFB is represented by a fi- 
nite number of box states. It has been demonstrated by 
explicit calculations for weakly bound nuclei [1, [iJl that 
such a box discretization is very accurate when compared 
to the exact results. 

In many practical applications involving complex ge- 
ometries of nucleonic densities in two or three dimen- 
sions, such as those appearing in nuclear fission and fu- 
sion or weakly bound and spatially-extended systems, it 
is crucial to consider large coordinate spaces. At the 
same time, employed lattice spacings should be suffi- 
ciently small to provide good resolution and numerical 
accuracy. As a result, the size of the discretized contin- 



uum space may often become intractable. This is also the 
case for calculations employing multiresolution pseudo- 
spectral methods [S^ which effectively invoke enormous 
continuum spaces. Therefore, it becomes exceedingly im- 
portant to develop methods allowing precise treatment 
of HFB resonances and non-resonant quasi-particle con- 
tinuum without resorting to the explicit computation of 
all box states. Our paper is devoted to this problem. 
Namely, we study the effect of high-energy, non-resonant 
quasi-particle continuum on HFB equations and observ- 
ables and propose an efficient scheme to account for those 
states. The technique, based on the local-density approx- 
imation to HFB-Popov equations for Bose gases [36|] , was 
previously used to approximate the continuum contribu- 
tion in solving the Bogoliubov de-Gennes equations for 
cold Fermi gases 37]. Here, we extend this method to 
the nuclear Skyrme HFB approach. 

We also study quasi-particle resonances and devise 
techniques to isolate them and estimate their widths us- 
ing C? discretization. Several approximate methods are 
examined. In atomic physics, the so-called stabilization 
method has been widely used to precisely calculate reso- 
nance widths p38l - l41| . A modified stabilization method, 
based on box solutions, was developed .40.1 a nd used to 
study single-particle resonances in nuclei [i^l. Here, we 
demonstrate that the stabilization method also works re- 
liably for quasi-particle resonances of HFB. Besides the 
stabilization method, a straightforward smoothing and 
fitting method that utilizes the density of box states is 
proposed and tested. Finally, we assess the quality of 
the perturbative expression 6] for deep-hole resonance 
widths. 

This paper is organized as follows. Section [H] briefly 
reviews the properties of HFB equations. In Sec. IIIII 
we apply the local-density approximation to account for 
the high-energy quasi-particle continuum. We propose 
and test a hybrid HFB strategy that makes it possible 
to solve HFB equations assuming a low energy cutoff, 
thus appreciably reducing the computational effort. Sec- 
tion IIVI studies deep-hole HFB resonances with several 
methods: the smoothing and fitting method, the box sta- 
bilization method, and the perturbation expression. The 
applicability of each technique is examined. Section |V] 
contains illustrative examples for weakly bound Ni iso- 
topes. Finally, conclusions are given in Sec. IVII 



II. COORDINATE-SPACE HFB FORMALISM 



and lower components of quasi-particle wave functions, 
respectively; Ei is the quasi-particle energy; and A is the 
Fermi energy (or chemical potential). For bound systems, 
A < and the self-consistent densities and fields are lo- 
calized in space. In our work, we shall limit discussion 
to local mean fields. To relate better to other papers, in 
the following we should use the notation A(r) = /i(r) for 
the pairing potential. 

For \Ei\ < —A, the eigenstates of ^ are discrete and 
Vi{r) and Ui(r) decay exponentially. The quasi-particle 
continuum corresponds to \Ek\ > —A. For those states, 
the upper component of the wave function always has a 
scattering asymptotic form. By applying the box bound- 
ary condition, the continuum becomes discretized and 
one obtains a finite number of continuum quasi-particles. 
In principle, the box solution representing the continuum 
can be close to the exact solution when a sufficiently big 
box and small mesh size are adopted. 

The solution of Eq. ([T]) in non-spherical boxes is a dif- 
ficult computational task. The recently developed par- 
allel 2D-HFB solvers utilizing the B-spline technique of- 
fer excellent accuracy when describing weakly bound nu- 
clei and large deformations [l8|, |32| . To solve the HFB 
equations in a 3D coordinate space is more complicated, 
but the development of a general-purpose 3D-HFB solver 
based on multiresolution analysis and wavelet expan- 
sion [32, Sa is underway. With the 2D-HFB box solver 
HFB- AX [33, the HFB equation is solved by discretiz- 
ing wave functions on a 2D lattice. The discretization 
precision depends on the order of B-splines, the maxi- 
mum mesh size, and the box size. Using this solver, one 
can obtain an extremely dense quasi-particle energy spec- 
trum by adopting a large box size. For example, there 
are about 7,000 states below 60 MeV if a square box of 
40x40 fm is used. In this case, the 2D box solution cor- 
responds to a high resolution of discretized continuum 
states. 



III. THOMAS-FERMI APPROXIMATION TO 
HIGH ENERGY HFB CONTINUUM 

In this section, the Thomas-Fermi (TF) approximation 
pBa . |37| is applied to study high-energy HFB continuum 
contributions. Calculations are carried out within the 
Skyrme HFB approach. The accuracy of this method is 
tested by comparing with the results obtained using the 
ftfll box discretization method. 



In this section, we briefiy recall general properties of 
the HFB eigenstates. The HFB equation in the coordi- 
nate space can be written as: 



Contribution from high-energy continuum 
states 



h 
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-h + X 



Uiir) 



= E,, 






(1) 



where h is the HF Hamiltonian; h is the pairing Hamil- 
tonian as defined in Ref. 71]; Ui and Vi are the upper 



The TF approximation to the HFB continuum was 
originally applied in the context of ID Fermi gases p37j . 
Following the method of [3a, |33 , we derive expressions 
for the continuum components to nucleonic densities due 
to high-energy scattering states within the Skyrme HFB 



approach. 

The high-energy HFB wave functions u{r) and v{r) 
can be approximated by |15l - [l7l l36l | : 

u(r) ^ u{p, r)e*'*'^('-) , w(r) -^ v(p, r)e''"'"*('-) , (2) 

where V0(r) = p. In this approximation, the deriva- 
tives of u{p, r) and v{p, r), and the second derivatives of 
(/)(r) are neglected, as well as the derivative terms of the 
effective mass and spin-orbit. The latter terms have neg- 
ligible effects on high-energy states, and they are us ually 
excluded in the pairing regularization procedure Il5l - ll7l | . 
The corresponding HFB equation ([T]) can be written as 



—— — h V(r) — A I uip, r) + A(r)v(p, r) 

2M* ^ ' J "-f^ ' \ ! \i-^ I ^g^^ 

= E{p,r)u{p,r), 



— — — -|- V{r) — A u(n, r) — A(r)u(p, r) 

2M* ^ "^ ' J y^' ' \ J \i^, J ^g^^ 

= ~E{p,r)v{p,r), 



where V{r) is the Skyrme HE potential; M*{r) is the 
density-dependent effective mass; and E{p,r) denotes 
the local quasi-particle energy. By introducing the HF 
energy enpiP, »") relative to the chemical potential. 



eHF(p,r) 



n2p2 



+ V{r)-X, 



2M*{r) 
E{p, r) takes the familiar BCS-like form: 



(4) 



Consequently, the contribution to the particle density 
from the high-energy continuum states is given by 



Pcir)=f^{l-'-f)e[E-E^], 



(8) 



where Ec is the quasi-particle energy cutoff above which 
the TF approximation is applied. In Eq. ([8|), the 
sum/integral symbol denotes the summation over the dis- 
cretized continuum box states or, alternatively, the inte- 
gration in the momentum space if the HFB equations 
are solved with the outgoing or scattering boundary con- 
ditions. Similarly, the high-energy continuum contribu- 
tions to the kinetic density, pairing density, and pairing 
potential are given respectively by: 



Tcir) = 

Pcir) = 
Ac(r) = 



^^'1-'J^)q[E-E^], 



VriUr) 



f p 

Vpair{r)pc{r) 
-Vpa^r{r)^^e[E^E,] 



Er 



(9) 
(10) 

(11) 



By separating the continuum contribution from the equa- 
tion ^ , we see that the TF procedure [Sa, 1231 is formally 
equivalent to the pairing regularization scheme with an 
effective pairing interaction Ve//(r) [l5l - [l7| : 



K//(t-) Vpair{r) 



Y,^Q[E-E^ 



(12) 



^(P,r-) = J4p(p,r)-|-A2(r). 



(5) 



In this work we consider the contact pairing interaction 
with the density-dependent pairing strength Vpair{i')- 
The resulting pairing potential A(r) can be written as: 



A(r) = Vpa^r{r)p{r) = \v^ 



1 



2po 



Kr), 



(6) 



where Vq is the pairing interaction strength; po is the 
saturation density 0.16 fm^'^; and p{r) and p{r) are the 
particle density and pairing density, respectively. The 
pairing potential © corresponds to the so-called mixed 
pairing interaction |43|]. 

The normalization condition in the p space, \u{p, t)P+ 
\v{p,r)\'^ = 1, implies that ^ 



z ^ 1 /, _ eHF(p,r) 



u; = -^ii 



E{p,r) 



gHF(P,r-) 
E{p, r) 



(7a) 



(7b) 



The expressions ([Sl llip can be written in a compact 
form by replacing the momentum sum/integral by an in- 
tegral over quasi-particle energy space between Ec and 
the maximum cutoff energy considered E„i : 



Pcir) = 

Tc{r) = 

Pc{r) = 
where 

piE,r) 



M*{r) 
27r2ft2 

M*{r) 
27r2ft2 



dE 



E 



dE 



M*{r)A{r) 



^y/E^-A^r) 

E 
^y/E^-A^r) 
^" p{E, r) dE 
E^ v/^2_A2(r) 




2M*{r) 
~l9 



^£;2_A2(r)-T/(r) + A . (14) 



The effective pairing strength (fT2|) becomes 



1 



1 



M*{r) 



Veffir) Vpa^r{r) 2Tl^fl^ J e^ y^ E^ - A^ (r) 



p{E,r)dE 



(15) 



To examine the quality of the TF approximation for 
high-energy continuum states, we compared it with the 




FIG. 1: (Color online) Continuum contributions in the energy 
range from Ec=40 MeV to 5^=60 MeV to (a) particle, (b) 
kinetic, and (c) pairing densities (see Eqs. (| 8110(1 ) of ™Zn 



obtained from the TF approximation and discretized box so- 
lutions (symbols) of hfb-AX with the SkM* energy density 
functional and mixed pairing interaction. See text for details. 



results of box discretization calculations obtained with 



for ™Zn with the SkM* 



the HFB solver hfb-AX j 

energy density functional (EDF) [4J] and mixed pair- 
ing interaction. The pairing strength is taken to be 
Vq = —234.85 MeV fm~'^ that is adjusted so as to repro- 
duce the neutron pairing gap of ^^"^Sn. The calculations 
have been carried out with the B-spline order of 12; the 



mesh size 0.6 fm, and the box size 24x24 fm. The nucleus 
^"^Zn is predicted to be spherical and has non-vanishing 
pairing in both protons and neutrons. 

Figure [1] displays the continuum contributions to par- 
ticle, kinetic, and pairing densities due to unbound states 
from Ec =40 MeV to E„i =60 MeV. (In this energy win- 
dow, the discretized HFB continuum contains no deep- 
hole resonances.) They were obtained from Eqs. P^ 
and from discretized hfb-AX solutions. It can be seen 
that the two methods produce very close continuum con- 
tributions to the local densities. For the neutrons, the 
continuum densities are mainly concentrated at the nu- 
clear surface. The proton densities have a more pro- 
nounced volume character. We found that this difference 
is mainly due to different pairing potentials A(r) and 
depends weakly on mean-field potentials V{r). The con- 
tinuum kinetic and pairing densities have similar shapes 
to the continuum particle densities. It can be seen, how- 
ever, that the continuum contributions to pairing densi- 
ties are larger than to continuum particle densities by two 
orders of magnitude. This is to be expected: the HFB 
continuum is strongly affected by the pairing channel [7| . 
Similar conclusions have been obtained in Ref. |20| . 

It is to be noted that the kinetic and pairing density 
integrals in Eq. (J13p are divergent and a finite upper limit 
for the integration Ec must be taken. In practice, the 
dependence on Ec can be avoided by adopting the pairing 
regularization pS]) . As discussed in Refs. (iSl - llTt , results 
obtained with regularized pairing are independent on the 
cutoff energy Ec and are very close to the results of a 
pairing renormalization procedure adopted in this work. 
Moreover, as also pointed out in Refs. [4^|4g' the kinetic 
energy density Tc has the same type of divergence as the 
pairing density pc, and the sum of kinetic and regularized 
pairing energies converges. 



B. The hybrid HFB strategy 

We have seen that the Thomas- Fermi approximation to 
the non-resonant continuum gives results very close to the 
box discretization approach. This suggests a hybrid HFB 
strategy to separately treat the deep-hole resonances 
and the high-energy non-resonant continuum. Figure [21 
schematically displays the HFB quasiparticle spectrum. 
As discussed earlier, the bound HFB solutions exist only 
in the energy region Ei ^ —A. The quasi-particle contin- 
uum with Ei > —A consists of non-resonant continuum 
and quasi-particle resonances. Among the latter ones, 
the deep-hole states play a distinct role. In the absence 
of pairing, a deep-hole excitation with energy Ei corre- 
sponds to an occupied HF state with energy Si = —Ei. 
If pairing is present, it generates a coupling of this state 
with unbound particle states with e^ « —Ei that gives 
rise to a quasi-particle resonance with a finite width. 

The low-energy continuum with —X<Ei^Ec con- 
sists of many resonances that have to be computed as 
precisely as possible. Therefore, to this region we apply 
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FIG. 2: (Color online) Schematic picture of the HFB quasi- 
particle spectrum. In the energy region Ei ^ —A, quasi- 
particle states are bound. In the hybrid HFB strategy, the 
quasi-particle continuum with Ei > —A is divided into the 
low-energy continuum with —X<Ei^Ec, which is treated 
by means of the box discretization, and the high-energy con- 
tinuum with Ec < Ei ^ Em that consists of non-resonant 
continuum treated by means of the TF approximation and 
several deep-hole states. See text for details. 



the box discretization technique. The high-energy con- 
tinuum with Ec < Ei ^ Em can be divided into the 
non-resonant part and several deep-hole states. While 
the non-resonant continuum can be integrated out by 
means of the TF approximation, deep-hole states have 
to be treated separately. Indeed, as discussed in Sec. |IV] 
below, deep-hole resonances are sufficiently narrow to be 
considered as a separate group of states. 

To this end, we diagonalize the HF Hamiltonian, 



(/^-A)^;r(0=^.«^W, 



(16) 



to obtain wave functions v^^{r) and single-particle en- 
ergies Ei of deep-hole states. The function vf^{r) is a 
very good approximation to the lower HFB component 
Vi{r). In the BCS approximation, assuming the state- 
dependent pairing gap (A^ 0,1131), 



A,;- A{r)[v^'{r)]'dr, 



the effective BCS occupation factor is: 



vl 



£k 



VM 



^l 



(17) 



(18) 



Consequently, v^ (r) is normalized to the BCS occu- 
pation vf. The corresponding quasi-particle energy is 
related to the single-particle energy through Ei = —£i. 

The hybrid HFB strategy is based on combining the 
box solutions below some low energy cutoff E^, the deep- 
hole solutions, and high-energy TF continuum. In this 



FIG. 3: The binding energy of ^"Zn as a function of the low- 
energy quasi-particle neutron cutoff Ec in the hybrid HFB 
method. The HFB-AX value corresponds to Ec=60MeV+\. 



way, the total particle density (as well as other HFB den- 
sities) of the even-even system can be split into three 
parts: 



Be 



p(r)= ^2K(r)p+ Y.2\v 

Ei>0 Eih>E, 



dh 



(r)p-Hp,(r), (19) 



where the continuum particle density pc of states between 
Ec and Em is given by Eq. (J13p . In practical applica- 
tions, the maximum cutoff energy Em is usually taken as 
60 MeV— A. The continuum contribution to the particle 
number, 



Nr. 



Pc{r)dr, 



(20) 



is always included to meet the particle number equation. 

We tested this hybrid HFB strategy to calculate the 
binding energy of ^"Zn at different cutoff values Ec- The 
HF equations (|16|) are solved using hfb-AX. Generally, 
the deep-hole single-particle energies Ei obtained from 
Eq. (|T6)) are very close to the hole-like solutions of the 
full HFB diagonalization. In ^"Zn, for example, the 
quasi-particle energy of the Isi/2 neutron is Ei=38.071 
MeV, while the corresponding HF single-particle energy 
is -£^=38.055 MeV. 

Figure 13] shows the total binding energy as a function 
of the neutron cutoff Ec- It is seen that the total en- 
ergy is perfectly stable for Ec + A„>30MeV, and it is 
equal to the binding energy obtained by means of the 
full bock discretization {Ec = Em)- At a low value of 
Ec + A„=15 MeV, the error of the hybrid method is only 
about 100 keV. At even lower values of cutoff, the TF 
approximation deteriorates rapidly [iTj . Note that the 
pairing strength in the hybrid HFB should not be renor- 
malized with Ec- The choice of the cutoff Ec is deter- 
mined by positions of deep-hole levels; this information 



can be obtained by solving the HF problem. The hybrid 
strategy can also be modified to work with the pairing 
regularization procedure [l3|- However, as seen in Fig.|3l 
this technique cannot be used with a very low cutoff Ec 
where resonances are densely populated. 



IV. PROPERTIES OF DEEP-HOLE 
RESONANCES OF HFB 

In this section, we study several approximate methods 
to calculate widths of HFB resonances. Calculations are 
performed for ^°Zn using the SkM* EDF and the mixed 
pairing interaction. 



A. The smoothing and fitting method 

From the box discretization of hfb-AX, one obtains a 
finite and very dense spectrum of continuum states. The 
deep-hole states are no longer isolated as in the BCS 
approximation; they become fragmented due to the pair- 
ing coupling with the neighboring particle- like continuum 
and acquire a decay width. The energy distribution of 
the occupation numbers 



\vi{r)f dr 



(21) 



has roughly the Breit-Wigner shape. Figure |4] displays 
occupation probabilities for the J7'^=l/2"'", l/2~, and 
3/2^ discretized neutron quasi-particle states in ^''Zn as 
a function of quasi-particle energy Ei. Note that the 
angular momentum projection f2 and parity tt are good 
quantum numbers in hfb-AX. It is apparent that the vf 
distributions of lsi/2i lP3/2i and lpi/2 deep-hole states 
have resonance-like structure. For spherical nuclei, de- 
generate resonances with different values oi\fl\ belonging 
to the same shell are expected to have the same width, 
and this is indeed the case for both magnetic substates of 
lp3/2. Breit-Wigner resonance-like structures in occupa- 
tion probabilities have also been predicted by the recent 
continuum HFB calculations [20|. In the Green function 
approach [^, [20], the occupation probability is related 
to the continuum level density, which corresponds to a 
Breit-Wigner shape for resonances [41|. In the present 
work, we propose a straightforward smoothing and fit- 
ting method to estimate resonance widths from discrete 
v'f distributions. 

To extract resonance parameters from the discrete dis- 
tribution of vf, we first smooth it using a Lorentzian 
shape function u;(.t) = 1/[2tt{x^ + 1/4)] (til,!!!]. The oc- 
cupation numbers vf of states above the Fermi level are 
smoothed out by means of folding with w{E/T), 
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FIG. 4: (Color online) Occupation numbers of the discretized 
neutron quasi-particle continuum states with Q^ — 1/2"*" 
(top), Q,'' = 1/2" (middle), and Q,'' = 3/2" (bottom) calcu- 
lated for ^"Zn with hfb-AX. Since ^"Zn is spherical, angular 
momenta of discretized resonance states are marked. The cor- 
responding Breit-Wigner envelopes are indicated by dashed 
lines. Note that the envelopes for IP3/2 magnetic substates 
are identical. 



HE) = E 



vf (E-E., 

r V r 



(22) 



where F is a smoothing parameter. 



The smoothed Breit-Wigner distribution is given by a 
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FIG. 5: (Color online) Smoothed occupation numbers (|22p 
and the corresponding smoothed Breit-Wigner distribution 
(gSJ for the Q." = 1/2+ neutron HFB resonances in ™Zn. 



compact expression [4l|, |48| : 



vUE)= ^{(E-Er? 



} 



(IT 1 ,2^;,^ r^ 



(23) 



-^(^-^0/3 In 



4^2 ^ p2 



4£;2 + r2 



where 

/3 = [(ii;-i?,)'+(r-r,) 2/4] "' [(i?-£;.)2+(r+ro V4] "'• 

(24) 
In the smoothing and fitting method, we fit the 
smoothed v'^{E) distribution ([22)1 using a two-term ex- 
pression: 



v^{E)^v''vl{E)+vl{E). 



(25) 



The first term is the smoothed Brcit-Wigner resonance 
shape. The background contribution v'^{E) is assumed to 
be a slowly changing Fermi-Dirac function ajieT^I^ -I- 1) 
characterized by parameters a and h. In order to deduce 
the resonance energy E^-, width F^, and the occupation 
number v^iv^ =^1), the method of least squares is used. 
The results of such a procedure for the Jl'^ = 1/2+ neu- 
tron HFB resonances in ™Zn are illustrated in Fig.O the 
fitted Breit-Wigner resonances agree very well with the 
smoothed occupation numbers. This demonstrates that 
the discretized occupation numbers have a Breit-Wigner 
shape. One has to note, however, that for threshold reso- 
nances (such as those close to E = —A), the distribution 
can strongly deviate from Breit-Wigner 1131] . 



As it has been pointed out in Ref. |4l|, the fitting 
curve could be dependent on the smoothing parameter 
F. We checked that such a dependence is not significant 
as far as the resonance width is concerned if the value 
of F around 0.8 MeV is taken. Generally, one should 
use F ^ F^. A very small smoothing parameter is not 
sufficient to smooth out the finite discretization effects. 
On the other hand, a large smoot hing parameter can un- 
derestimate the resonance width [41[. Numerical errors 
can also arise from the fitting procedure if several res- 
onances are overlapping. The discretization can yield 
a crude representation of occupation distribution if too 
small a box is used [20]. A high discretization resolu- 
tion is of particular importance when it comes to narrow 
resonances. Another advantage of the direct smoothing- 
fitting method discussed here is that it can be applied to 
resonances in deformed nuclei. 



B. The box stabilization method 

The stabilization method ^] is an C? method, and 
has been used to obtain precisely the resonance energy 
Ej. and widths F^. in atomic [3^, Ho, S^j and nuclear 
|50l ISlj physics. Based on the box solutions, the HFB 
resonances are expected to be localized solutions with 
energies weakly affected by changes of the box size. The 
stabilization method allows to obtain the resonance pa- 
rameters from the box-size dependence of quasi-particle 
eigenvalues. 

To this end, one introduces the continuum level density 
^{E\ 



g{E) - 1:t[5{E -H)- 5{E - Ho)] 



(26) 



where Hq is the non-interacting Hamiltonian of the sys- 
tem. The continuum level density is related to the phase 
shift by 



S{E) =TT g{E')dE'. 
Jo 



(27) 



In a modified stabilization method, one can obtain 
the phase shift using the box-discretized continuum [40| ■ 
Firstly, we compute quasi-particle energies Ej{L) with 
differentbox size L using the spherical HFB solver hf- 
BRAD [30] . Figure M shows that the discretized quasi- 
particle energies generally smoothly decrease as the box 
size increases. Around the resonances, however, the ener- 
gies are fairly constant. Starting from a sufficiently large 
box Lq, calculations are done on a grid with the spacing 
SL: 



L = Lo + Ai, AL = {M - l)6L. 



(28) 



In this way, eigenvalues Ej are stored in arrays Ej (L). By 
using the Akima interpolation |52l [. we obtain the values 
of Lj{E) corresponding to a box size with which the j'^ 
eigenvalue equals to E. 




FIG. 6: (Color online) The neutron discretized quasi-particle 
£—0 spectrum of ^"Zn as a function of the box size L. Calcu- 
lations have been carried out with the spherical HFB solver 
HFBRAD [33]. Unlike the non-resonant box continuum, reso- 
nances are practically L-independent. 



In the next step, the phase shift is obtained from [40l |: 

S{E) = nN{E) + ^ ^(Lo + AL - L,{E)), (29) 
j 

where N{E) is the number of eigenvalues for vifhich 
Ej{Lq) < E. To obtain the resonance energy E^ and 
the width F^, we fit the phase shift S{E) by using the 
expression: 

5{E) = arctan /'^(^__^^ + g^^^^), (30) 



where 6b{E) represents the smoothly changing back- 
ground contribution parametrized as: 



Sb{E) =a + b{E-Er). 



(31) 



Figure [7] shows the phase shifts of the neutron lsi/2 
resonance obtained with different step numbers M. It 
can be seen that the background becomes more smooth 
as the step numbers increase. Generally, a large AL is 
necessary to smooth out the background, so as to reduce 
the fitting errors in a and b in Eq. (j3ip . In our calculation, 
we also find that for very narrow resonances, a small 
value of 6L is important for the interpolation precision. 
In principle, as 



SL -^ and AL — ;■ 00, 



(32) 



one can get very accurate resonance parameters. In prac- 
tice, however, calculations with a large box and small step 
sizes are very expensive even with a fast hfbrad solver. 
The eigenvalue spectrum of hfbrad is sparse, resulting 
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FIG. 7: The phase shifts 5{E) of the neutron lsi/2 state of 
^'^Zn obtained with the stabilization method with different 
iteration steps (|28[) . Here, SL=0.06 fm and Lo—24: fm. 



in statistical errors in the calculations of phase shift. For 
low-energy HFB resonances, this situation is even worse 
since these eigenvalues change very slowly as the box size 
increases, as shown in Fig. [HI In this case, the calculated 
background is not sufficiently smooth. Currently, the 
stabilization method is limited to spherical boxes as in 
HFBRAD. The extension of the stabilization method to 
deformed cases will be the subject of future work. 



C. The perturbative expression 

Assuming that the pairing coupling can be treated per- 
turbatively, the width of a deep-hole HFB resonance can 
be given by the Fermi Golden rule 0, Ifl, 01 ■ 



F = 2^|(uo£|A|i- 



HP I 



(33) 



where A is the pairing potential; v^^ is the bound wave 
function of the HF potential corresponding to the single- 
particle energy —E+X, and uqe is the scattering solution 
of the HF equation with energy ii^ -I- A. 

At a large distance, the scattering wave function uqe 
has the usual asymptotic form: 



1 

UOE = - 

r 



2m 
hPnk 



(cos diFi{Tj,kr)+ sin diGi{r],kr)) , (34) 



where Fi and Gi are the regular and irregular Coulomb 
functions, respectively; di is the phase shift; and k — 
^j2m{E + \)/h?. The scattering wave function is found 
with the help of Ref . [53 . In the calculation of scatter- 
ing wave functions, h^ /2m is taken as a constant 20.734 
MeV, to avoid problems related to the density-dependent 
effective mass. 



TABLE I: Energies (in MeV) and widths (in keV) of deep- hole HFB resonances in ^"Zn calculated with the box stabilization 
method (box), smoothing- fitting method (smf), and perturbative expression (per). See text for more details. 







Neutron 






Proton 




states 


E 


i box A smf 


J- per 


E 


A box t smf 


J- per 


lSl/2 


38.075 


1.23 0.80 


0.98 


35.147 


0.64 0.54 


0.40 


lP3/2 


28.801 


3.42 4.1 


0.04 


25.424 


3.22 2.04 


2.27 


lPl/2 


26.043 


6.01 7.0 


0.19 


22.851 


2.04 1.57 


1.26 


1^5/2 


18.987 


13.1 10.6 


22.7 


15.320 


0.12 


0.09 



Table U displays the widths of high-energy deep-hole 
states in ™Zn calculated with the smoothing-fitting 
method, box stabilization method, and perturbative ex- 
pression. It is seen that the widths obtained with the 
smoothing-fitting method and box stabilization method 
are always fairly close. For protons, all three methods 
yield similar results. However, the perturbative expres- 
sion predicts widths that are too small for the IP3/2 and 
lPi/2 neutron states. The difference is due to a surface- 
peaked neutron pairing potential that is very different 
from that of protons (see FigH]). Indeed, for surface- 
peaked pairing, the overlap between uqe (which exhibits 
rapid spatial oscillations), w^^, and A(r) is too sensitive 
to small changes in uqe- For more discussion on limi- 
tations of the perturbative expression (p3| . see also the 
schematic model analysis in Ref. [6|. 

Compared to other deep-hole resonances in Table U 
the width of the proton ld5/2 state is very small. This 
is because the upper component wave function of this 
i = 2 state is quasi bound, due to the confining effect 
of the Coulomb-plus-centrifugal barrier (see Fig. [51). For 
such narrow resonances, the stabilization method cannot 
easily be applied. 



APPLICATIONS TO WEAKLY BOUND 
NUCLEI 



The effects of continuum are expected to increase when 
moving towards the neutron drip line |54l |. Especially 
for very weakly bound nuclei exhibiting halo structures, 
there exists a strong interplay between pairing and con- 
tinuum |55l-l57j . Here, we investigate the role of contin- 
uum contributions in the neutron-rich Ni isotopes. These 
nuclei have been studied by the Gamow HFB (GHFB) 
method and the exact quasi-particle resonance widths of 
^•^Ni are available [8]. The GHFB work provides an ex- 
cellent benchmark for our approximate resonance widths. 

As in Ref. 8J, the calculations have been carried out 
with the SLy4 EDF [53| and the surface pairing interac- 
tion f43|, with the strength Vb=-519.9 MeV fm'^. In 
Figini we show the high-energy continuum contributions 
to the neutron pairing densities in 84,86,88,90]\^j^ obtained 
using the TF approximation p^ . It can be seen that 
the continuum pairing density acts in the surface region, 
and increases as nuclei become less bound. This is consis- 
tent with the earlier findings [7^] that continuum becomes 
important for nuclei close to the drip line. 
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FIG. 8: (Color online) A special kind of deep-hole HFB proton 
resonance, in which the upper component of the HFB wave 
function is quasi bound, due to the confining effect of the 
Coulomb-plus-centrifugal barrier . 



FIG. 9: (Color online) The continuum contributions to the 
neutron pairing densities of neutron-rich 84,86,88,90-|s^j isotopes, 

=60 MeV 



calculated with Eq. ^3^ with Ec=30 MeV and E^-- 
and mixed pairing interaction. See text for details. 
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FIG. 10: (Color online) Similar to Fig.|9]but with the volume, 
surface, and mixed pairing interaction in **Ni. 



Figure [TU] illustrates the interaction dependence of the 
pairing-continuum effect. Namely, it displays the neu- 
tron continuum pairing densities for ^*Ni with the vol- 
ume pairing (Vb=— 185.026 MeV fm~^), surface pairing 
as in Figini and mixed pairing (Vb=— 284.36 MeV fm~^). 
(The pairing strengths Vq have all been adjusted to re- 
produce the neutron pairing gap of ^^"^Sn.) It is seen 
that the high-energy continuum contribution to pairing 
density strongly depends on character of pairing interac- 
tion [6, 7]. In particular, in the case of surface pairing, 
the continuum contribution is remarkably larger than for 
other two pairing functionals, indicating its very different 
behavior in weakly bound nuclei [12| . Actually, as it has 
been discussed 32] the surface pairing is essential for the 
existence of bound ^"Ni. 

Table |ll] displays the widths of neutron resonances in 
^•^Ni calculated using different techniques. In the cal- 
culations with the box stabilization method, we used 
Lo=22 fm, M=400, and SL=OM fm, see Eq. ^ . 
Figure [11] illustrates the quality of calculations for the 
neutron lp3/2 resonance. It is seen that the fitted phase 
shift agrees well with that obtained from the stabilization 
method. Systematically, the box stabilization method 
predicts slightly larger widths as compared to GHFB. 
This is consistent with findings of Ref. [39| where the 
stabilization method generally overestimates the widths 
by 10%. In particular, it is shown that the widths of 
very narrow resonances are largely overestimated. The 
width of the lsi/2 state is so narrow that it is beyond 
the applicability of the stabilization method. The very 
narrow lg9/2 state belongs to the class of special HFB 
resonances of FiglS) Among the resonances, the lg7/2 
and lhii/2 states are also single-particle resonances in 
the Hartree-Fock theory Q. Other than that. Table HIl 
demonstrates that the stabilization method works well 
for all the HFB resonances except for extremely narrow 
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FIG. 11: (Color online) The phase shift of the neutron IP3/2 
resonance in ^°Ni using Eqs. (|3T]| and (|29)l. 



Within the smoothing-fitting method, the quasi- 
particle energy spectrum is obtained by hfb-AX by tak- 
ing a large box of 38x38 fm. The widths given by the 
smoothing-fitting method agree with the exact numbers 
within a factor of two. We have found that the fitting 
precision is compromised when several resonances over- 
lap. For the low-energy resonances, the total occupation 
numbers v'^ in Eq. (j25p are very small and can induce 
additional fitting errors. Besides, as we have discussed 
earlier, it is not proper to fit the occupation probabilities 
of low-energy resonances near the Fermi energy using a 
Breit-Wigner shape. In spite of all those reservations, the 
precision of the smoothing-fitting method is quite satis- 
factory. 

For neutron HFB resonance widths in Ni isotopes, 
our calculations predict that generally the widths would 
slowly increase as the drip line is approached. For exam- 
ple, the widths of the neutron IP3/2 state in ^^Ni, ^^Ni, 
and ^°Ni calculated by the stabilization method are, re- 
spectively, 25.97 keV, 28.25 keV, and 30.84 keV. This 
is consistent with Fig. |9j the widths grow with the in- 
creased pairing-continuum coupling. 



VI. CONCLUSIONS 

In this work, we performed a comprehensive study of 
quasi-particle continuum within the HFB theory. The 
purpose of this investigation is twofold. Firstly, we tested 
a truncation scheme based on the Thomas-Fermi approx- 
imation to limit the continuum space in realistic calcu- 
lations carried out in huge configuration spaces (or large 
spatial boxes) that yield huge amounts of discretized un- 
bound states. Secondly, we studied properties of HFB 
resonances, including deep-hole states. We compare sev- 
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TABLE II: Energies (in MeV) and widths (in keV) of HFB neutron resonances in ^''Ni, calculated using the box stabilization 
(box) and smoothing-fitting (smf) methods. They are compared to the GHFB solutions of Ref. [g]. 



states 


E, 


Tbox 


Tsmf 


Tghfe 


lSi/2 


51.419 


- 


l.le-3 


1.09e-3 


lP3/2 


40.588 


30.84 


20.17 


27.28 


lPl/2 


38.770 


32.26 


34.67 


27.14 


14/2 


29.039 


1.31 


1.37 


0.78 


1^3/2 


25.017 


25.44 


23.08 


22.57 


2Si/2 


24.319 


50.36 


40.87 


46.00 


1/7/2 


17.554 


401.79 


413.04 


397.37 


2P3/2 


12.538 


499.98 


471.22 


490.56 


1/5/2 


10.981 


672.44 


651.12 


645.64 


2pi/2 


10.816 


440.76 


376.19 


404.30 


lff9/2 


6.519 


1.56 


0.52 


0.81 


2^5/2 


3.270 


221.08 


105.17 


194.18 


2^3/2 


2.310 


611.50 


643.88 


560.61 


lff7/2 


3.348 


69.09 


75.13 


63.61 


I/11I/2 


5.527 


162.46 


173.66 


131.78 



eral methods to estimate resonance width and discuss 
their strengths and weaknesses. 

The TF approximation to the high-energy contin- 
uum states in the hybrid HFB variant, in which deep- 
hole states and non-resonant continuum are separately 
treated, has been found very effective for it fully repro- 
duces results obtained with the full box discretization. 
The high-energy non-resonant continuum has been found 
to have a similar spatial behavior for particle, kinetic, and 
pairing densities. This distribution is mainly determined 
by the pairing potential. The continuum contribution to 
the pairing density is substantial for weakly bound nu- 
clei and it has appreciable spatial extension. The hybrid 
method will be useful in reducing the computational cost 
of 3D coordinate-space HFB calculations. 

The HFB quasi-particle resonance is unique to the 
HFB theory and is a fascinating topic in its own right. 
We examined three approximate methods to study the 
resonance widths based on HFB box solutions. In the 
smoothing-fitting method, resonance parameters are ob- 
tained by fitting the smoothed occupation numbers ob- 
tained from discretized solutions. The box stabilization 
method is based on the fact that quasi-particle ener- 
gies of continuum states change with the box size. By 
comparing with the exact Gamow HFB results obtained 
by imposing outgoing boundary conditions [8|, we have 
demonstrated that the stabilization method works fairly 
well for all HFB resonances, except for the very narrow 
ones. The smoothing-fitting method is also very effec- 
tive and can easily be extended to deformed cases. The 
perturbative Fermi golden rule is found to be unreliable 
for calculating widths of neutron resonances. The only 
exceptions are narrow metastable deep-hole states such 
as high-£ states and low- lying proton resonances. 

Illustrative examples have been provided for the drip- 
line Ni isotopes. We found that continuum densities 
strongly depend on the density dependence of pairing 
interaction. In particular, surface pairing produces very 



large continuum pairing densities. The obtained neutron 
widths of ^°Ni are generally larger than that of stable nu- 
clei. The presence of broad quasi-particle resonances in 
weakly bound nuclei suggests that quasi-particle contin- 
uum plays an important role in the description of excited 
states. In addition, we expect that the determination of 
neutron resonance widths can be useful to estimate neu- 
tron emission half-lives of excited states above the neu- 
tron emission threshold. 

In summary, we have demonstrated that one can im- 
plement powerful approximations to incorporate the vast 
quasi-particle continuum space without explicitly impos- 
ing scattering or outgoing boundary conditions. We 
expect our work can also be useful in the context of 
continuum- QRPA applications. The approximate tech- 
niques used in this study have demonstrated the precision 
of the box discretization in representing the continuum 
and deep-hole resonances, especially important for nuclei 
near drip lines where continuum effects are large. 
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